clear all
close all
clc

addpath('BC')
addpath('ex32_toolbox_local')
addpath('../exercise23/EFG')
addpath('../exercise23/EFGChino')

RHOx=2.6;%%%2.6;
PENALTY=10^4  /8%80000;


%% EFG options (Machine precision) this is for compatibility!
optEFG.dim   = 2;
optEFG.verb  = 0;
optEFG.grad  = 1;           % Computation of the Gradient 0:OFF 1:ON
optEFG.hess  = 0;           % Computation of the Hessian  0:OFF 1:ON
optEFG.TolNR = 1.e-10;
optEFG.knn   = 0;
optEFG.Tol0  = 1.e-06;
optEFG.gamma = 1.8;         % value of gamma to compute LME

% Nodal Integration options
% cubature order: 1 (1 GPts), 2 (3 GPts), 4 (6 GPts), 5 (7 GPts), 6 (12 GPts)
optGL.orderGL = 30;

%% Matrial Parameters
L  = 4;         % length
D  = 1;         % height
P  = -1;     % external load
E  = 10^03;     % Young modulus
nu = 0.3;       % Poisson coefficient
parameters.L  = L;
parameters.D  = D;
parameters.E  = E;
parameters.nu = nu;
parameters.P  = P;
Lx = L;
Ly = D;

%% Convervence curves
m      = 4;
m_id   = (2:m);      %index of the L2 error norm to perform the fit
h_n    = zeros(m,1);
L2_err = zeros(m,1);
%NChen  = [5 9 17 33];
NChen  = [5 9 17 33 33];

for n=1:m
  fprintf(1,'n= %d  init',n);
  % Node Points
  Ny = ceil(NChen(n)/2);
  Nx = NChen(n)*4-3;

  h  = L / (Nx - 1);
  optEFG.rho = RHOx * h;
  optEFG.spacing = h;
  

  % This line is important
  center  = 0.5*[Lx 0];
  x_nodes = UniformGrid2D(Nx, Ny, Lx, Ly, center);

  nPts  = length(x_nodes);
  h_n(n)= 1/sqrt(2*NChen(n)*Nx);
  
%   ids   = (1:nPts);
%   id_bd = [ids(x_nodes(:,1)==0) ids(x_nodes(:,2)==0) ...
%     ids(x_nodes(:,1)==Lx) ids(x_nodes(:,2)==Ly)];
%   id_bd = unique(id_bd);

  % Sample points and gauss weights
  [x_samples w_samples] = MakeGLSamples2D(x_nodes, optGL);
    
  %% COMPATIBILITY - Nodes thermalization, nodes and samples adjacency structures are computed
  %[beta_n range_n] = NodalThermalization(x_nodes, optEFG);
  optEFG.beta   = optEFG.gamma/(h*h);
%  optEFG.beta_n = beta_n;
%  optEFG.range_n= range_n;
  

  %% for adjacency calculation, using four times the rho
  range_n = ones(size(x_nodes,1))*4*optEFG.rho;
  optEFG.range_n= range_n;
  
  %% The basis functions are computed
  % adjacency structure with the nearest neighbors nodes to each sample point
  % nodal shape parameter
  s_near = SamplesAdjacency_efg(x_nodes,x_samples,range_n);
  % Local-max entropy basis functions computation
  optEFG.s_near = s_near;
  outEFG = wrapperCHINO(x_nodes, x_samples, optEFG);
  p_samp  = outEFG.p_samp;
  dp_samp = outEFG.dp_samp;

  PENALTY = PENALTY / h
  %% The displacement field is computed
  [u_h condK] = CantileverBeam_Solver_BC_penalty_COMPLETE_BAR(dp_samp,s_near,x_nodes,x_samples,w_samples,parameters,optEFG, PENALTY);

  %% L2 norm
  u_num  = reshape(u_h,2,length(x_nodes))';
  u_num  = SamplesSolution(u_num,s_near,p_samp);
  u_sol  = ex32_AnalyticalSolution(x_samples,P,D,L,E,nu);
  L2_err(n) = sqrt(sum(sum((u_num-u_sol).^2,2).*w_samples))/sqrt(sum(sum(u_sol.^2,2).*w_samples))
end

%% Analysis of Results
fit_num  = polyfit(log10(h_n(m_id)),log10(L2_err(m_id)),1);
disp('- - - - - - - - - - - - - - -');
fprintf ( 1, 'EFG  : m=%7.4g  b=%7.4g\n', fit_num(1),  fit_num(2));
% Plots
h_line=[0.01 0.1]*2;
err_line=[0.0001 0.01]*2;
fit_line=polyfit(log10(h_line),log10(err_line),1);
fprintf ( 1, 'line : m=%7.4g  b=%7.4g\n', fit_line(1),fit_line(2));
figure(1);clf
loglog(h_n,L2_err,'ro-',h_line, err_line,'k-', ...
  'LineWidth',2,'MarkerSize',6)
xlabel('h','fontsize',14,'fontweight','b')
ylabel('Relative Error','fontsize',14,'fontweight','b')
legend('EFG','Location','Best','Orientation','horizontal')
title(strcat('Relative Error in Norm L_2 -- \gamma_{EFG}=',num2str(optEFG.gamma)),'fontsize',16,'fontweight','b')
set(gca,'XTick',0.02:0.02:0.1)
%xlim([0.001 0.11])


% Node and Sample Points
figure(2)
plot(x_nodes(:,1),x_nodes(:,2),'or',...
  x_samples(:,1),x_samples(:,2),'.b');
axis equal
legend('Node points','Sample points');


% Deformed and undeformed (numerical and analytical) configurations
figure(3)
u_sol = ex32_AnalyticalSolution(x_nodes,P,D,L,E,nu);
u_sol = u_sol+x_nodes;
u_num = reshape(u_h,2,length(x_nodes))';
u_num = u_num+x_nodes;
plot(x_nodes(:,1),x_nodes(:,2),'.r')
axis equal
hold on
plot(u_num(:,1),u_num(:,2),'xb')
hold on
plot(u_sol(:,1),u_sol(:,2),'.g')
hold on
legend('Undeformed','Deformed (numerical)','Deformed (analytical)');



%% Calcula los valores del gradiente de la soluciona proximada
% calcula los valores de u en todos los puntos de gauss
u_num  = reshape(u_h, 2, length(x_nodes))';
u_num  = SamplesSolution(u_num, s_near, p_samp);
% quiero mostrar: para cada sample point xs, el valor de 
gu = zeros(2,2,size(x_samples,1));
for i = 1 : size(x_samples,1);
    nears = s_near{i};
    gradN = dp_samp{i};
    for b = 1 : size(nears);
        gu(:,:,i) = gu(:,:,i) + tensorProduct(gradN(b,:), u_num(i,:), 2);
    end;
end;

k=1;
for j = 1 : size(x_samples,1);
    if (abs( x_samples(j, 2) - x_samples(1, 2) ) < 1.e-7);
        indXG(k) = j;
        k=k+1;
    end;    
end;

gUxx  = gu(1,1, indXG);
gUxy  = gu(1,2, indXG);
gUyx  = gu(2,1, indXG);
gUyy  = gu(2,2, indXG);
gUxx  = reshape(gUxx, numel(indXG),1);
gUxy  = reshape(gUxy, numel(indXG),1);
gUyx  = reshape(gUyx, numel(indXG),1);
gUyy  = reshape(gUyy, numel(indXG),1);


